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FAST  ERROR-FREE  ALGORITHMS 
FOR  POLYNOMIAL  MATRIX  COMPUTATIONS 


John  S.  Baras,  David  C.  MacEnany  and  Robert  L.  Munach 
Systems  Research  Center 
University  of  Maryland 
College  Park,  MD  20742 

Abstract 

Polynomial  matrices  over  fields  and  rings  provide  a  unifying  framework  for  many 
control  system  design  problems.  These  include  dynamic  compensator  design,  infinite  di¬ 
mensional  systems,  controllers  for  nonlinear  systems,  and  even  controllers  for  discrete  event 
systems.  An  important  obstacle  for  utilizing  these  powerful  mathematical  tools  in  practical 
applications  has  been  the  non-availability  of  efficient  and  fast  algoritms  to  carry  through 
the  precise  error-free  computations  required  by  these  algebraic  methods.  Recently,  with 
the  advent  of  computer  algebra  this  has  become  possible.  In  this  paper  we  develop  highly 
efficient,  error-free  algorithms,  for  most  of  the  important  computations  needed  in  linear 
systems  over  fields  or  rings.  We  show  that  the  structure  of  the  underlying  rings  and 
modules  is  critical  in  designing  such  algorithms.  We  also  discuss  the  importance  of  such 
algorithms  for  controller  synthesis. 


1.  Introduction 

The  theory  of  polynomial  matrices  [8,21,23]  plays  a  key  role  in  the  frequency-domain 
approach  to  the  synthesis  of  multi-input/multi-output  control  and  communication  systems 
[13,24,25].  Examples  include  coprime  factorizations  of  transfer  function  matrices,  canon¬ 
ical  realizations  obtained  from  matrix  fraction  descriptions,  design  of  feedback  -compen¬ 
sators  and  convolutional  coders,  and  the  analysis  of  quantization  effects  in  linear  systems. 
Typically,  such  problems  abstract  in  a  natural  way  to  the  need  to  solve  systems  of  gen¬ 
eralized  Diophantine  equations,  e.g.,  the  so-called  Bezout  equation  [6,15,19,22].  These 
and  other  problems  involving  polynomial  matrices  require  efficient  polynomial  matrix  tri- 
angularization  procedures  [16],  a  result  which  is  not  surprising  given  the  importance  of 


1 


matrix  triangularization  techniques  in  numerical  linear  algebra.  There,  matrices  with  en¬ 
tries  from  a  field  can  be  triangularized  using  some  form  of  Gaussian  elimination.  However, 
polynomial  matrices  have  entries  from  a  polynomial  ring,  an  algebraic  object  for  which 
Gaussian  elimination  is  not  defined.  For  matrices  with  entries  from  a  polynomial  ring 
which  is  Euclidean — the  kind  encountered  most  often  in  control  theory  applications — 
triangularization  is  accomplished  instead  by  what  is  naturally  referred  to  as  Euclidean 
elimination.  Unfortunately,  the  numerical  stability  and  sensitivity  issues  of  Euclidean 
elimination  are  not  well  understood  and  in  practice  floating-point  arithmetic  has  yielded 
poor  results.  At  present  r  a  reliable  numerical  algorithm  for  the  triangularization  of  poly¬ 
nomial  matrices  does  not  exist. 

This  paper  presents  algorithms  for  polynomial  matrix  triangularization  which  entirely 
circumvent  the  numerical  sensitivity  issues  of  floating-point  methods  through  the  use  of 
exact,  symbolic  methods  from  computer  algebra  [5,14,20].  Often  one  encounters  the  com¬ 
ment  that  since  in  practical  problems  the  numerical  coefficients  are  not  known  precisely, 
computer  algebra  methods  are  at  a  disadvantage  from  a  practical  point  of  view.  However 
this  is  a  misconception.  Whether  we  know  the  coefficients  accurately  or  not  is  not  the 
issue.  The  real  issue  in  using  these  algebraic  methods  is  to  what  extend  we  can  perform 
the  required  computations  within  the  accuracy  of  the  model  data.  Existing  floating  point 
methods  are  poor,  highly  sensitive  and  often  lead  to  large  errors;  essentially  since  they 
suffer  from  the  same  problems  as  computing  zeroes  of  polynomials.  The  use  of  exact, 
error-free  algorithms  guarantees  that  all  calculations  are  accurate  to  within  the  precision 
of  the  model  data— -the  best  that  can  be  achieved.  Furthermore  one  can  calculate  with  such 
algorithms  the  exact  sensitivities  involved  and  therefore  judge  appropriately  the  confidence 
one  should  place  on  the  results.  Previous  computer  algebra  algorithms  for  polynomial  ma¬ 
trix  problems  appearing  in  control  systems  have  been  reported  in  [11].  Their  performnce 
was  very  slow  even  on  small  size  problems.  We  place  emphasis  on  efficient  algorithms  to 
compute  exact  Hermite  forms  of  polynomial  matrices.  The  triangular,  or  more  correctly, 
trapezoidal  Hermite  form  is  defined  for  any  matrix  with  entries  from  a  principal  ideal  ring 
[21,23].  Such  matrices  arise  in  many  practical  problems  in  communications  and  control. 
Here  we  shall  focus  on  matrices  having  entries  which  are  polynomials  with  rational  coef- 
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ficients.  An  important  aspect  of  the  exact  triangularization  of  such  matrices  involves  the 
choice  of  arithmetic.  We  consider  the  tradeoffs  between  rational  and  integer  arithmetic 
and  choose  the  latter.  This  choice  leads  us  to  consider  algorithms  for  the  division  of  poly¬ 
nomials  over  a  unique  factorization  domain  (UFD).  The  standard  algorithm  for  this  task 
is  well-known  [3,4,7,18]  and  defined  more  generally  for  polynomials  with  coefficients  from 
any  commutative  ring  with  identity.  This  algorithm  is  well-suited  to  the  scalar  problem  of 
GCD  computation  of  polynomials  over  UFDs  since  it  avoids  the  computation  of  GCDs  of 
the  coefficients.  In  the  context  of  polynomial  matrix  triangularization  however,  it  becomes 
unavoidable  to  exploit  the  richer  structure  of  the  coefficient  ring:  the  fact  that  GCDs  are 
defined  on  a  UFD.  As  a  result  we  present  an  alternative  to  the  standard  algorithm  special-, 
ized  to  polynomials  over  UFDs  but  enjoying  a  certain  optimality  property  which  is  crucial 
to  the  efficiency  of  matrix  triangularization  procedures. 

We  have  implemented  algorithms  to  compute  exact  Hermite  forms  of  polynomial  ma¬ 
trices  in  the  MACSYMA  and  Mathematica  computer  algebra  languages.  We  have  also 
written  a  suite  of  auxiliary  programs  which  call  on  these  triangularization  procedures  in 
order  to  perform  the  more  high-level  tasks  arising  in  the  frequency- domain  approach  to 
control  system  synthesis.  We  conducted  simulations  with  MACSYMA  code  running  on 
Texas  Instruments  Explorer  II  and  give  performance  results  for  the  triangularization  of 
polynomial  matrices.  The  algorithms  presented  here  have  time  performance  which  is  three 
orders  of  magnitude  better  than  algorithms  known  earlier  to  control  engineers. 

2.  Facts  And  Terminology  of  Polynomials  and  Polynomial  Matrices 

In  this  section  we  draw  upon  standard  material  from  modern  algebra  [10,12],  Denote 
by  Q[s]  the  ring  of  polynomials  in  the  indeterminate  ‘s’  with  coefficients  drawn  from 
the  field  of  rational  numbers,  Q.  The  subring  Z[s]  of  Q[s]  results  when  the  polynomial 
coefficients  are  restricted  to  lie  in  Z ,  the  ring  of  integers.  The  leading  coefficient  of  a 
polynomial  is  the  nonzero  coefficient  of  its  highest  degree  term.  By  convention,  the  leading 
coefficient  of  the  zero  polynomial  is  defined  to  be  one,  and  its  degree  is  taken  to  be  —  oo. 
If  the  leading  coefficient  of  a  polynomial  is  one,  then  the  polynomial  is  said  to  be  monic. 
If  a(s)  in  Z[s]  is  a  polynomial  of  degree  zero,  then  a{s)  is  said  to  be  a  unit  of  Q[s],  i.e., 
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(a(s))-1  is  in  Q[s].  Obviously,  the  units  of  Q[s]  are  precisely  the  units  of  Q,  namely, 
Q  \  {0}.  A  polynomial  a(s)  in  Z[s\  is  called  primitive  if  its  coefficients  are  relatively  prim,e 
in  Z ,  i.e.,  if  the  greatest  common  divisor  of  its  coefficients  is  a  unit  of  Z,  namely  ±1. 
For  any  a(s )  in  Z[s],  there  exists  a  non-zero  scalar  ca  in  Z,  unique  up  to  its  sign,  and  a 
primitive  polynomial  pa(s )  in  Z[s],  such  that  a(s)  =  ca  -pa(s).  The  pair  (ca,pa(s))  is  called 
a  content-primitive  factorization  of  a(s);  with  slight  imprecision  ca  is  called  the  content  of 
a(s )  and  pa(s )  its  primitive  (with  respect  to  ca ).  A  collection  of  polynomials  in  Z[s ]  having 
contents  which  are  relatively  prime  is  said  to  be  relatively  primitive.  A  polynomial  p(s) 
divides  a  polynomial  q(sj-,  written  p(s)|^(s),  if  there  exists  c(s)  such  that  p(s)  c(s)  =  q(s). 
A  common  divisor ,  CD,  c(s)  of  {p;(s)  :  i  =  1, . . . ,  n}  is  a  polynomial  such  that  c(s)|{pt(.s)}. 
A  greatest  common  divisor ,  GCD,  g(s)  of  {p,(s)}  is  a  CD  of  {p;(s)}  such  that  c(s)|^(s) 
for  any  other  CD,  c(s),  of  {p^(s)}. 

Denote  by  M[Q[s]]  the  collection  ofmxn  matrices  with  entries  from  Q[s].  Similarly,. 
M[Z[s]]  will  denote  the  subset  of  M[Q[s]]  when  the  entries  are  resticted  to  lie  in  Z[s\.  We 
call  A(s)  in  M[Q[s]]  a  polynomial  matrix.  Letting  R  denote  either  Z[s]  or  Q[s],  A(s)  in 
M[R]  is  said  to  be  nonsingular  if  A($)  is  square  (m  —  n)  and  det  A(s)  is  not  the  zero 
polynomial.  A  nonsingular  polynomial  matrix,  U(s)  in  M[R)  having  a  determinant  which 
is  a  unit  of  R  is  a  unit  of  M[R\  and  is  said  to  be  unimodular  (with  respect  to  M[R ]). 
In  this  case  therefore,  U(s)~x  —  Adj  U(s)/  det  U(s)  is  itself  a  polynomial  matrix  and  also 
unimodular  with  respect  to  M[R].  When  the  ring  R  also  happens  to  be  a  field,  the  concepts 
of  unimodularity  and  nonsingularity  conincide,  but  keep  in  mind  that  in  the  general  ring 
case  it  is  unimodularity  that  coincides  with  the  usual  notion  of  invertibility. 

A  row  of  a  polynomial  matrix  A(s)  in  M[Z[s]]  is  said  to  be  primitive  if  its  polynomial 
entries  are  relatively  primitive;  A(s)  is  said  to  be  row  (left)  primitive ,  if  every  row  is 
primitive.  For  any  A(s)  in  M[Z[s]],  there  exists  a  diagonal  matrix  Ca  in  M[Z]  and  a 
row  primitive  matrix  Pa(s)  in  M[Z[,s]]  such  that  A(s)  —  Ca  •  Pa{s)-  This  we  call  a  left 
content-primitive  factorization  of  A(s).  The  diagonal  elements  of  Ca  are  called  the  row 
contents  of  the  respective  rows  of  A(s).  By  analogy  with  the  scalar  case,  a  content-primitive 
factorization  is  obviously  unique  only  up  to  the  choice  of  the  signs  of  the  row  contents. 

Two  polynomial  matrices  Ai(s)  and  A^fs)  in  M[Q[s]]  (or  polynomial  vectors  if  m  =  1 
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or  n  =  1)  are  said  to  be  linearly  dependent  if  there  exists  a  coefficient  of  proportionality 
a(s)/b(s )  in  the  field  of  rational  functions  such  that  -4i(s)  =  a(s)/b(s)A2(s).  Hence  Ai(s) 
and  A2(s)  are  linearly  dependent  if  there  exist  a(s)  and  6(s)  in  Q[s]  such  that  b(s)Ai(s)  — 
a(s)d.2(s)  =  0.  This  leads  to  the  simple  generalization  that  a  finite  collection  of  polynomial 
matrices  is  linearly  dependent  if  and  only  if  there  exists  a  nontrivial  linear  combination 
of  them  employing  only  polynomial  coefficients  which  equals  the  zero  polynomial  matrix. 
With  this  in  mind,  elementary  row  and  column  operations  for  polynomial  matrices  in 
M  [Q[s]]  take  the  form:  ' 

1.  interchange  of  any  two  rows; 

2.  addition  to  any  row  a  polynomial  multiple  of  any  other  row; 

3.  multiplication  of  a  row  by  a  nonzero  rational. 

Performing  any  of  the  above  operations  on  an  identity  matrix  results  in  an  elementary 
matrix  E(s)  in  M[Q[s]].  Clearly,  each  such  E(s )  is  unimodular  as  is  its  inverse  and  also 
the  product  of  any  number  of  elementary  matrices.  Conversely,  every  unimodular  matrix  is 
a  product  of  elementary  matrices.  Two  polynomial  matrices  >1(5)  and  B(s )  are  said  to  be 
row  equivalent  if  each  can  be  obtained  from  the  other  using  a  finite  sequence  of  elementary 
row  operations.  In  other  words,  they  are  related  by  left  multiplication  with  a  unimodular 
matrix  U(s)  such  that  B(s)  =  U(s)  A(s)  and  A(s)  —  U~1(s)  B(s). 

Every  m  x  n  polynomial  matrix  A(s)  in  M[Q[s]]  is  row  equivalent  to  an  upper  tri¬ 
angular  form  (or  upper  trapezoidal  form  if  m  ^  n).  Therefore,  it  can  be  reduced  by  a 
sequence  of  elementary  row  operations  to  an  upper  triangular  (trapezoidal)  matrix  Ha(s) 
in  M[Q[s}].  In  fact,  there  exists  a  unimodular  matrix  U(s )  such  that  U(s )  A(s)  =  Ha(s ) 
with  Ha(s)  satisfying  the  following  conditions: 

1.  Each  entry  below  the  diagonal  is  identically  zero; 

2.  Each  nonzero  diagonal  entry  has  degree  greater  than  the  entries  above  it; 

3.  Each  diagonal  entry  is  monic. 

We  say  that  Ha(s)  is  a  column  monic-Hermite  form  of  A(s).  To  state  conditions  more 
precise  than  these  concerning  the  structure  of  Ha(s )  requires  the  notion  of  the  rank  of 
a  polynomial  matrix  which  we  omit.  This  omission  is  forgivea.ble  in  that  rank  should 
not  enter  into  an  exact  triangularization  algorithm.  Indeed,  exact  rank  information  for 
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polynomial  matrices  is  often  obtained  most  efficiently  through  triangulariza,tion. 

A  column  integral- Hermite  form  can  be  defined  in  terms  of  the  column  monic-Hermite 
form.  Letting  Ha(s)  denote  a  column  monic-Hermite  form  for  A(s)  in  Af  [Q[.s]],  multiply 
each  row  of  H/ i(s)  with  the  respectively  smallest  positive  integer  such  that  the  matrix 
H'a(s )  so  obtained  is  in  Af[Z[.s]].  Clearly,  H'A(s )  is  row  primitive  and  row  equivalent  to 
A(s).  Conversely,  suppose  that  one  is  given  H'A(s )  satisfying  conditions  (1)  and  (2)  above 
which  is  row  primitive  and  row  equivalent  to  A(s).  Divide  each  row  of  H'A(s)  by  the  leading 
coefficient  of  the  polynomial  on  the  diagonal  of  the  respective  row  and  call  the  matrix  so 
obtained  Ha(s).  Then  clearly  there  existsT7(s)  unimodular  such  that  U(s)A(s )  =  Ha(s) 
and  Ha(s)  is  a  monic-Hermite  form  of  A(.s).  This  concept  of  column  integral-Hermite  form 
gives  a  triangular  form  in  M[Z[s]]  for  each  matrix  in  M[Q[s]].  If  A(s)  is  nonsingular  then 
it  can  be  shown  that  its  monic-Hermite  form  is  unique  and  therefore  its  integral-Hermite 
form  is  also  unique. 


3.  Triangularizing  Polynomial  Matrices 

The  upper  triangularization  of  matrices  with  entries  from  a  field  using  a  sequence  of 
non-singular  elementary  row  operations  plays  a  key  role  in  the  application  of  the  theory 
of  vector  spaces.  Likewise,  the  upper  triangularization  of  matrices  with  entries  from  a 
ring  using  a  sequence  of  unim,odular  elementary  row  operations  plays  a  key  role  in  the 
application  of  the  theory  of  vector  modules.  Recall  that  the  axioms  defining  a  vector 
module  (or  matrix  m.odule  or  simply  module)  over  a  ring  R  are  the  same  as  those  defining 
a  linear  vector  space  except  that  linear  combinations  employ  scalar  coefficients  from  a 
general  ring  R  without  the  restriction  that  R  also  be  a  field.  Thus  a  vector  space  over  a 
field  F,  V(F),  is  precisely  the  matrix  module  over  the  ring  F ,  Ad(F).  We  only  deal  with 
msitrix  modules  over  commutative  rings  (with  identity)  so  there  is  no  need  to  distinguish 
module  handedness.  We  employ  the  notation  M[A](R)  to  denote  the  m  x  n  matrix  module 
over  the  ring  R  with  matrix  entries  from  the  additive  abelian  group  A;  for  the  case  M[R](R) 
we  have  been  writing  M[R\.  Using  this  notation  we  point  out  that  M[Q](Z )  is  a  module 
but  that  M[Z](Q )  is  not.  We  also  warn  the  reader  that  AL[Z[s]]  is  not  a  submodule  of 
M[Q[s]\:  a  module  AL[R](/?)  is  said  to  be  a  submodule  of  a  module  M[A](R)  if  B  is  an 
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additive  subgroup  of  A  and  for  each  b  in  B ,  rb  is  in  B  for  all  r  in  R.  Thus,  M[^[s]]  is  a 
submodule  of  M[Q[s]](Z[s]). 

Computing  triangular  (trapezoidal)  forms  of  matrices  can  be  accomplished  on  any  ma¬ 
trix  module  of  the  form  M[R\  where  R  is  an  integral  domain ,  i.e.,  a  commutative  ring  with 
identity  having  no  zero  divisors  [1,2,17].  However,  this  cannot  in  general  be  accomplished 
using  only  unimodular  (i.e.,  invertible)  operations.  Nevertheless,  the  transformation  to  an 
upper  triangular  form  using  unimodular  elementary  row  operations  can  be  performed  quite 
straightforwardly — in  theory  at  least— on  any  matrix  with  entries  from  a  type  of  integral 
domain  called  a  Euclidean  ring,  for  instance  on  a  ma-trix  from  M[Q[s]].  The  key  feature 
that  Euclidean  rings  enjoy  is  the  Euclidean  division  property  which  we  state  for  Q[s].  Given 
polynomials  a(s),  h(s)  in  Q[s]  there  exist  two  unique  polynomials,  the  quotient  q(s)  and 
the  remainder  r(s),  such  that  b(s)  =  q(s)  a(s)  -(-  r(s)  and  degr(s)  <  dega(s).  Often  this 
property  is  stated  under  the  assumption  that  dega(s)  <  deg  b(s)  since  the  proof  for  the 
case  deg  b(s)  <  dega(s)  is  trivial:  choose  r(s)  =  a(s)  and  q(s )  =  0.  The  Euclidean  division 
lemma  for  a  general  Euclidean  ring  R  has  a  statement  which  differs  from  the  above  only  in 
that  polynomial  degree,  the  ring  valuation  for  Q[$],  is  replaced  by  the  ring  valuation  for  R] 
for  R  =  Z  (the  case  considered  by  Euclid)  the  ring  valuation  is  integer  absolute  value.  The 
fact  that  the  inequality  on  the  degrees  of  a(s)  and  r(s)  is  strict  allows  one  to  introduce  a 
zero  into  a  polynomial  matrix  using  elementary  operations.  As  an  example,  consider  an 
attempt  to  introduce  a  zero  into  the  (2, 1)  position  of  a  2  x  2  polynomial  matrix.  Suppose 
one  has  found  q(s)  and  r(s)  as  in  the  division  lemma  such  that  6(s)  =  </(s)  a(s )  +  r(s)  and 
then  computes, 


(  1  °\(a(s)  c(s)\  =  (a(s)  <s)  \ 

V~?(5)  vUs)  dls)J  W5)  d(s)  -  q(s)c(s)  J 

using  an  obviously  unimodular  pre-multiplication.  Notice  that  the  degree  of  the  (2,1) 
entry  has  been  decreased  by  at  least  one.  By  interchanging  the  rows,  the  same  procedure 
can  now  be  repeated  on  the  resulting  matrix  and  iterated,  each  step  reducing  the  degree 
of  the  (-2, 1)  entry  in  view  of  the  strict  degree  inequality.  This  process  consists  entirely  of 
unimodular  premultiplications  and  clearly  cannot  continue  indefinitely  without  yielding  a 
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constant  in  the  (2,1)  entry.  If  the  constant  is  not  zero,  one  final  row  exchange  and  the 
obvious  elementary  row  operation  will  introduce  a  zero  into  the  (2, 1)  position.  Moreover, 
one  can  show  that  the  resulting  polynomial  in  the  (1,1)  position  will  be  the  GCD  of 
the  original  polynomials  a(s)  and  6(s).  The  use  of  this  process  to  introduce  zeroes  into 
polynomial  matrices  we  call  Euclidean  elimination  by  analogy  with  Gaussian  elimination. 


4.  Integer  vs  Rational  Arithmetic 


The  above  example  also  serves  to  illustrate  the  fact  that  rational  arithmetic  Is  costly. 
For  instance,  to  calculate  the  coefficients, of  polynomials  of  the  form  d(s)  —  q(s)  c(s )  with 
c,  <i,  q  in  Q[s],  one  encounters  the  generic  computation  a  +  /3  y  with  oi, /?,  7  in  Q.  If  these 
rationals  are  expressed  as  ratios  of  integers  a  —  jd  =  7  =  all  reduced  to 

lowest  terms,  then 


a  +  7  8  = 


Na  D?  +  N0  N't  Da 
Da  DP  D~i 


This  computation  requires  six  integer  multiplications,  one  integer  addition  and  the  calcu¬ 
lation  of  a  GCD.  Although  there  are  more  efficient  methods  [17],  it  remains  a  fact  that 
rational  arithmetic  is  computationally  expensive,  due  in  large  part  to  the  need  for  GCD 
calculations.  On  the  other  hand,  if  it  can  be  arranged  so  that  a,  (3  and  7  are  all  integers, 
then  the  same  computation  obviously  requires  only  two  integer  multiplications,  one  integer 
addition  and  no  GCD  calculation.  Thus,  our  goal  is  to  carry  out  matrix  triangularization 
on  M[Q[s]]  using  only  integer  arithmetic.  Clearly,  by  multiplying  each  row  of  any  A(.s)  in 
M[<5[s]]  by  a  large  enough  integer,  the  denominators  of  every  coefficient  of  every  entry  of 
A(s)  can  be  cancelled  and  such  a  diagonal  operation  is  certainly  unimodular  in  M[Q[s]]. 
Again,  this  computation  can  be  arranged  more  efficiently  but  because  it  involves  a  fixed 
overhead,  assume  for  convenience  that  A(s)  is  given  in  M[Z[s]]. 

Unfort xmately,  this  creates  new  difficulties  because  Euclidean  elimination  is  not  de¬ 
fined  for  M[Z[s]]  since  Z[s\  is  not  a  Euclidean  ring.  For  instance,  it  is  easy  to  see  that 
the  remainder  of  two  polynomials  in  Q[.s]  with  integer  coefficients  has,  in  general,  rational 
coefficients;  consider  the  remainder  of  2s  after  division  by  3s  —  1.  In  other  words.  Euclidean 
division  is  not  defined  for  Z[s\.  However,  Z[s\  is  an  instance  of  a  polynomial  ring  with 
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coefficients  from  a  commutative  ring  with  identity  and  for  such  a  ring  one  has  the  pseudo- 
division  lemma ,  a  natural  generalization  of  the  Euclidean  division  lemma.  Let  C  denote  a 
commutative  ring  with  identity.  Given  a(s)  and  6(s)  in  C[s]  with  deg  a(s)  <  deg  b(s)  there 
exist  two  polynomials,  the  pseudo -quotient  q(s)  and  the  pseudo-remainder  r(s),  such  that 
Lb(s)  =  q(s)  a(s)  +  r(s)  and  deg  r(s)  <  deg a(s)  where  the  premultiplier  L  =  %(gb  dega+1 
with  ao  denoting  the  leading  coefficient  of  a(s).  The  pseudo-quotient  and  pseudo-remainder 
are  unique  if  C  is  also  an  integral  domain.  The  proof  of  the  pseudo-division  lemma  yields 
a  division  procedure  called  pseudo-division  which  like  Euclidean  division  enjoys  the  all- 
important  strict  degree  reduction  property;, see  [18]  for  the  standard  pseu do- division  algo¬ 
rithm. 

Let’s  consider  an  example  in  which  we  wish  to  pseudo-divide  b{ s )  by  a(s)  where, 
b(s)  =  s8  +  /  -  3s4  -  3s3  +  8s2  +  2s  -  5 

and 

a(s)  =  3s6  +  5s4  -  4s2  -  9s  +  21. 

Applying  the  standard  pseudo-division  algorithm  one  obtains, 

27  6(a)  =  (9s2  -  6)  a(s )  +  (-15s4  +  3s2  -  9), 

i.e.,  L  =  38_6+1  =  27,  q(s)  =  9s2  —  6  and  r(s)  =  —15s4  +  3s2  —  9.  This  example  appears 
in  [18]  as  one  step  in  the  task  of  computing  the  GCD  of  b(s )  and  a(s).  The  next  step  is 
to  divide  out  the  content  of  r(s)  and  then  compute  the  GCD  of  a(s)  and  pr(s)  exploiting 
the  fact  that  gcd(6(s),  a(s))  =  gcd(a(s),pr(s)).  The  purpose  of  this  content  removal  is  to 
keep  the  size  of  the  coefficients  small  for  purposes  of  efficiency  in  succeeding  calculations. 
However,  consider  the  above  computation  in  the  context  of  a  matrix  triangularization — our 
2  x  2  example  will  suffice: 

(l  0\  / a(s )  c(s)\  =  / «(s)  c(s) 

\-q(s)  Lj\b(s )  d{s)J  \  r(s)  L  d(s)  -  q{s)c(s) 

In  this  situation  we  see  that  we  are  not  at  liberty  to  blindly  divide  the  entire  second  row 
by  the  content  of  r(s )  (or  any  integer  for  that  matter)  because  it  may  introduce  rational 
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coefficients  in  the  (2,  2)  entry  and  thereby  ruin  our  attempt  to  maintain  integer  arithmetic. 
However,  note  that  another  solution  to  the  above  pseudo-division  example  is, 

9  b(s )  =  (3s2  —  2)  a(s)  +  (—5s4  +  s2  —  3), 

i.e.,  L  =  27  is  not  necessarily  the  smallest  premultiplier  for  which  a  “pseudo-quotient” 
and  “pseudo-remainder”  exist.  Obviously,  in  the  matrix  case,  “L  =  9”  yields  better  results 
than  L  —  21  since  it  yields  smaller  coefficients  in  the  second  row.  Of  course  in  this  example 
the  difference  is  negligible,  however,  if  the  size  of  the  leading  coefficient  of  a(s)  is  large, 
the  difference  in  computational  burden  can  be  quite  substantial.  Moreover,  as  we  shall  see 
below,  keeping  the  size  of  all  coefficients  as  small  as  possible  is  a  primary  goal. 

5.  Pseudo-division  for  Polynomials  over  a  UFD 

It  is  apparent  that  there  are  smaller  (and  larger)  premultipliers,  L,  than  the  one 
defined  in  the  pseudo-division  lemma.  Now  the  pseudo-division  lemma  is  the  best  that  one 
can  do  in  general  for  polynomials  over  a  commutative  ring  with  identity.  But  be  aware  that 
the  concept  of  ‘smaller’  referred  to  in  the  pseudo-division  example  is  inherited  from  the  fact 
that  Z  is  also  a  unique  factorization  domain  (UFD).  Recall,  a  UFD  is  a  commutative  ring 
with  identity  which  has  no  zero  divisors  and  admits  prime  factorizations.  Let  U  denote 
a  UFD.  One  can  think  of  u  in  U  as  being  “smaller”  than  u'  in  U  if  u  is  a  divisor  of  u' . 
For  the  problem  of  pseudo-division  of  polynomials  a(s),  6(s)  in  U[s],  what  we  seek  is  the 
smallest  premultiplier  L*  in  U  such  that  if  there  exist  L  in  U  and  q,  r  in  U[s]  satisfying, 

L  b(s)  =  q(s)  a(s)  +  r(s)  and  degr(s)  <  dega(s), 
then  L*  divides  L  and  q*(s),  ?-*(s)  in  C7[s]  exist  such  that, 

L*  b(s)  =  q*(s)a(s)  +  r*(s)  and  degr*(s)  <  deg  a(s). 

A  representation  for  L*  is  easy  to  obtain.  Let  U'  denote  the  field  of  quotients  of  U 
and  consider  U  as  embedded  in  U1 .  The  field  of  quotients  of  a  ring  is  defined  for  any  ring 
which  is  an  integral  domain  and  thus  is  defined  for  any  UFD;  for  example  Q  is  the  field 
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of  quotients  of  Z.  Consider  «(.s),  b(s)  in  U[s]  and  therefore  in  U' [s]  in  the  sense  of  the 
embedding.  It  is  a  fact  that  U' [s]  is  a  Euclidean  ring  since  U'  is  a  field  and  polynomial 
rings  over  fields  are  Euclidean.  Therefore  by  the  Euclidean  division  lemma  there  exist  q(s), 
p(s)  in  I7'[s]  such  that, 

b(s )  =  g(s)  a(s)  +  p(s)  and  degp(s)  <  dega(s). 

Now  write  p(s)  =  n(s)/D  for  n(s)  in  U[s]  and  D  in  U  with  D  and  the  content  of  n(s) 
relatively  prime  in  U .  Then, 

D  b(s)  =  n(s)  a(s)  +  D  p(s)  and  deg  D  p(s)  <  deg  a(s). 

Clearly,  D  p(s)  is  in  U[s]  so  that  one  must  conclude  that  A*,  q*(s),  and  r*(s)  repectively 
differ  from  D,  n(s)  and  D  p(s)  by  at  most  a  unit  of  U.  Thus,  we  may  choose  L*  =  D 
and  this  suggests  one  method  to  compute  T*  and  then  by  uniqueness  q*(s)  =  n(s)  and 
r*(s)  =  D  p(s).  Letting  L  =  6-dcg  “+1  (assuming  deg  a(s)  <  deg  b(s))  and  making  an 
appeal  to  the  pseudo-division  lemma  it  is  not  hard  to  show  that  L*  =  L/  gcd(L,  cq)  where 
cq  denotes  the  content  of  the  pseudo- quotient  for  the  premultiplier  L.  Finally  note  that, 

ged  (L,cq)  =  gcd(gcd(a0,5o),gcd(a0,9i),. . . ,  gcd(a0,  qn-m))- 

Although  this  gives  us  a  representation  for  L*,  it  does  not  lead  to  the  most  efficient 
computation  for  L*  due  in  large  part  to  the  fact  that  the  coefficients  of  q(s)  are  already 
too  “big”.  Thus  what  we  want  is  an  efficient  algorithm  to  compute  L*,  q*(s)  and  r*(s). 
The  algorithm  given  next  computes  the  L*,  q*(s)  and  r*(s)  as  discussed  above.  It  is  a 
distinct  improvement  over  the  pseudo-division  lemma  given  in  [18]  in  that  it  computes 
with  smaller  numbers.  It  does  so  by  exploiting  the  richer  structure  of  polynomial  rings 
with  coefficients  from  a  UFD  but  at  the  cost  of  both  generality  and  GCD  calculations. 
However,  in  the  matrix  problems  we  consider  this  cost  is  unavoidable. 

Algorithm  M  -  Pseudo-division  of  Polynomials  over  a  UFD 

Given  two  nonzero  polynomials  b(s)  —  bosn  +  &isn_1  +•••  +  &«  and  a(s)  —  aosm  + 
cii-s7"'”1  4-  .  . .  +  am  in  U[s]  with  m  <  n,  this  algorithm  computes  the  smallest  A*,  pseudoi- 
quotient  q*(s ),  and  pseudo-remainder  r*(.s)  as  discussed  above.  It  computes  L*,  q*(s),  and 
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r*(s)  directly  by  computing  GCD’s  on  the  fly.  This  involves  “smaller”  numbers  than  first 
using  the  pseudo-division  algorithm  and  then  computing  GCDs  as  outlined  above.  Bigger 
numbers  cost  more  in  GCD  calculations  and  given  the  size  of  the  integers  encountered  in 
polynomial  matrix  computations,  e.g.,  easily  greater  than  1000  digits,  this  algorithm  can 
save  a  substantial  amount  of  time.  For  simpler  notation  we  drop  the  ‘asterisk’  subscript 
in  the  algorithm’s  definition. 

BEGIN: 

mnm  <—  min  (m,  n  —  m)  - 

g  < —  GCD(bo,ao) 
l  <-  a0/g 

L  4-  l 

b0  <-  b0/g 

For  i  —  1  thru  n  —  m  Do 

For  j  =  i  thru  n-mDo 
bj  < —  bj  *  l 
EndDo 

For  j  =  1  thru  min  (mnm,  n  —  m  —  i  +  1)  Do 
bj+i— l  *  bj+i— i  ~~  dj  *  bi 
EndDo 

g  GCD(b{,  a0) 
l  *—  ao /g 
L  4—  L  *  d 
bi  *-  bi/g 
EndDo 
END 

The  algorithm  terminates  with  the  first  n  —  m  +  1  coefficients  of  b(s)  overwritten  accord¬ 
ing  to  {£>o,  b\, . . . ,  £>„-m}  {<7o,  9i,  •  •  • , qn-m}  and  the  remaining  coefficients  over  written 

according  to  {bn_m+1, . . . ,  bn}  {r0, . . . ,  rm_x }. 

Algorithm  M  -  Proof  of  Correctness 
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The  informal  language  description  of  Algorithm  M  basically  implements  the  following 
recursion  for  k  =  0, 1, . . . ,  n  —  m, 

&<_1)(s)  =  !>M; 

St  =  gcd(a0,4t_1))i 

\ k  —  «0  /  9k  •) 

Pk\= 

b^k\s)^lkb{k~1\s)-pka(s). 

Observe  that  deg  b^k\s)  <  deg  for  k  =  0, 1, . . .  ,  n  —  m  because  lk  b^  ^  =  pk  a0 

(where  b^  of  course  denotes  the  leading  coefficient  of  b^).  Hence,  degb^n~m\s)  < 
dega(s).  From  the  algorithm’s  definition  we  see  that  r*(s)  =  b('n~m\s)  and  X*  =  ^Jfc- 

Solving  the  recursion  above  we  obtain 

b ^  ^('<’)  F*6(s)  —  (p0/l  ‘  ‘  ‘  ^n  —  ms  “b  ’  ’  ’  Pn—m  —  l^n  —  ms  T  Pn  —  ni)  a(s')- 

The  algorithm  therefore  yields, 


n—m 


<z*(5)  =  S 

fc=o 


PkS 


•> 


whence, 

X*  b(s)  =  q*(s)  a(s)  +  r*(s)  and  deg  r*(s)  <  deg  a(s), 

with  X*  in  U  and  q*(s),  r*(s)  in  £7[s].  Thus  the  algorithm  indeed  computes  a  valid  solution; 
next  we  show  that  it  is  optimal.  Suppose  there  exist  another  L  in  U  and  q(s),  r(&)  in  U[s] 
such  that, 

L  6(s)  =  qys)  a(s)  +  r(s)  and  deg  r(s)  <  deg  a(s). 

Then  by  commutativity  X*  L  b(s)  =  L  X*  b(s)  implies, 

(X</*(s)  —  X*</(s))a(s)  =  X*  ?'(s)  —  Xr*(s). 
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Since  there  are  no  divisors  of  zero  in  a  UFD  this  gives, 


deg(L  q*(s)  -  L *  q(s ))  +  deg  a(s)  =  deg(L*  r(s)  -  Lr*(s)). 

Since  deg(L*  r(s)  —  Lr*(s))  <  max{deg  r(s),  deg  r*(s)}  <  dega(s)  it  must  be  true  that, 


deg(Ag*(.s)  -  L*q(s))  =  -oo, 


and  therefore  X*  q(s)  —  Lq*(s).  By  equating  coefficients  we  obtain, 


.nr-o*.- 


Pk=L*qk  k  =  0,1,...  ,n  -m, 


so  that, 


<7*  ]^[  h  =  Pk  L  k  -  0,1,..., 


n  —  m. 


i=0 


For  k  =  0  we  get  /q <7o  =  Lpo  and  therefore  lo\Lpo-  However,  from  the  definition  of  the 
algorithm  Z0  and  p0  are  relatively  prime  in  U,  or  coprime ,  and  so  in  fact  Z0  j L.  For  k  —  1 
we  get  lohqi  =  Lpi  and  therefore  /i[(^-)pi.  Again,  by  construction  Zj  and  pi  are  coprime 
and  so  h\j^-  In  general  we  have  lk  and  pk  coprime  and  lkqk  —  ( lo.Llk  |  )Pk  so  that  for 
k  --  n  —  m  we  obtain, 

L 


In— m  I 


z0  •  •  •  z 


As  a  result  Zq  •  •  •  Zn_m  |L  and  therefore  L*|L. 


QED 


6.  Pseudo-Euclidean  Elimination 

The  introducton  of  a  zero  below  the  diagonal  of  a  matrix  A(s)  in  M[Z[s]]  can  now  be 
performed  using  Algorithm  M.  This  procedure  we  shall  call  pseudo- Euclidean  elimination 
(PSEE)  for  obvious  reasons.  Consider  triangularizing  the  matrix: 


1  s  s 

.■  1  ( )  —  |  4o5  — 105  10  3s2  “I-  s  10 

7  —  5s  6s2  —  1  4s2  —  10 


14 


Pseudo-Euclidean  elimination  yields, 


-1351755s3  +  1373940s2  -  5102550s  -  7152750 
-43605s3  +  77103s2  -  234189s  -  35190 

3706425s4  -  5202000s3  +  18532125s2  +  15671025s  +  7152750 

:  -  -  / 

and  this  illustrates  the  main  disadvantage  of  triangularization  on  M[Q[s]]  performed  over 
M \Z{s]\ — the  coefficient  growth  of  the  polynomials.  As  the  number  of  rows  and  columns 
in  the  matrix  increases,  this  coefficient  growth  continues  unabated  and  begins  to  erode 
the  advantage  of  using  integer  arithmetic.  One  approach  to  handle  this  new  source  of 
coefficient  growth  is  to  remove  the  content  of  the  current  row  after  each  pseudo-Euclidean 
division  step.  It  is  better  to  remove  the  row  content  as  soon  as  possible  in  this  way 
rather  than  waiting  due  to  the  cost  of  computing  GCDs  of  large  integers,  nevertheless,  we 
illustrate  row  content  removal  for  the  current  example.  Factoring  the  above  matrix  into  a 
left  content-primitive  form  Ca  H'a{s)  yields, 


7577325  0 

0  89145 

0  0 
V 


/  765  0  0  \  /  9905  0  -1767s3  +  1796s2  -  6670s  -  9350  \ 

0  9  0  0  9905  -4845s3  +  8567s2  -  26021s  -  3910  . 

\  0  0  65025/  \  0  0  57s4  -  80s3  +  285s2  +  241s  +  110  / 

The  superfluous  left  content  of  this  matrix  can  then  be  discarded  since  this  is  equivalent  to 
multiplying  A(s)  by  CA  1  thereby  keeping  the  coeffcient  size  to  a  minimum.  We  emphasize 
that  C a  is  unimodular  with  respect  to  M[Q[s]]  but  not  with  respect  to  M[Z[s]].  We  stress 
that  up  to  the  signs  of  the  entries  across  the  rows  H'A(s)  is  the  same  matrix  which  would 
have  resulted  had  we  employed  row  content  removal  after  each  pseudo-Euclidean  division 
step  and  that  this  is  the  more  efficient  strategy.  Note  that  the  above  polynomial  matrix 
H'a(s )  is  nonsingular  and  in  column  integral-Hermite  form  and  that  therefore  the  unique 
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column  monic-Hermite  form  of  A(s)  is  obtained  directly  from  H'A(s)  as, 

1  n  1767a3  ,  1796a2  1334a  1870 

U  9905  9905  1981  1981 

f|  1  969a3  ,  8567a2  _  26021a  _  782 

U  1  1981  '  9905  9905  1981 

0  0  _  80a!  +  5s2  +  241a  +  m 

We  see  that  PSEE  provides  an  efficient  triangularization  procedure  for  M[Zjs]]  but, 
strictly  speaking,  PSEE  modified  with  content  factorization  is  not  a  valid  triangularization 
procedure  for  M[Z[s]]  because  content  removal  is  not  a  unimodular  operation  in  M[Z[s] ]. 
On  the  other  hand,  augmenting  PSEE  with  content  factorization  is  a  unimodular  opera¬ 
tion  for  M[Q[s]]  and  yields  an  efficient  triangularization  procedure  for  M[Q[s]]  by  avoiding 
rational  arithmetic  while  maintaining  integers  of  the  smallest  possible  magnitude  through¬ 
out  an  elimination.  An  analogous  remark  can  be  made  about  any  pair  of  modules  of  the 
form  M[f7[s]]  and  Af[{7,[s]]  as  discussed  above. 
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7.  Algorithms  for  Triangularizing  Polynomial  Matrices 

Algorithm  T  -  Column- Oriented  Triangularization  of  Polynomial  Matrices 

Given  an  N  x  N  nonsingular  matrix  A  £  M[Z[s]],  this  algorithm  overwrites  A  with  a 
triangular  form  obtained  by  a  sequence  of  unimodular,  elementary  row  operations.  It 
avoids  rational  arithmetic  by  using  pseudo-division  as  defined  in  Algorithm  M  in  order 
to  achieve  maximum  computational  efficiency  with  minimum  coefficient  growth.  In  ad¬ 
dition,  it  further  inhibits  coefficient  growth  by  factoring  out  the  feow  content  after  each 
pseudo-Euclidean  division  step.  This  algorithm  operates  in  a  column  oriented  fashion  by 
successively  zeroing  out  the  entries  in  each  column  below  the  diagonal.  This-is  shown 
pictorially  below. 
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Assume  there  exists  a  pre-defined  function 

MinimumDegreeIndex(A ,  k )  :=  argmin{deg  A\^,  ■  ■  ■ ,  deg  Ajv,fc} 

which  returns  the  index  of  the  row  of  A  whose  kth  entry  is  a  non-zero  polynomial  of  lowest 
degree  among  the  rows  {k,  k  +  1,...,JV}.  If  Aktk{s)  —  ^k+ i,k(s)  —  -^N,k(s)  =  0,  then  it 
returns  —  oo,  the  degree  of  the  zero  polynomial  by  our  convention. 

BEGIN: 

For  k  =  1  thru  N-l  Do 

index  <—  MinimumDegreelndex^A,  k) 

If  index  ^  —oo  Then 

Akt.  Aindet,.  (exchange,  rows  k  and  index ) 

For  n  =  /,:  -|-  1  thru  N  Do  (zero  out  all  entries  in  column  k  below  Aktk) 
EndlessLoop 

num  <—  pseudo  —  -quotient(An  k , 
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denom  <—  pseudo  —  remainder^An^,  Ak,k) 

An<,  «—  denom  *  A„y  —  num  *  Ak,. 

An>,  <—  An,./GCD{conte?it(An,1), . . . ,  content{An^)} 

If  AUjk  =  0  then  exit  EndlessLoop 

An,.  Ak,. 

End  EndlessLoop 

EndDo  - 

End  If 
EndDo 
End 

It  is  emphasized  that  this  is  a  ‘paper’  algorithm.  In  fact,  the  working  code  based  on  the 
above  is  more  efficient  and  can  handle  singular  (we  have  not  defined  the  notion  of  rank 
of  a  polynomial  matrix  here)  and  non-square  matrices  and  the  entries  can  initially  belong 
to  Q[s],  However,  these  programming  considerations  just  complicate  matters  and  obscure 
the  basic  operation. 


Algorithm  P  -  Principal  Minor- Oriented  Triangularization  of  Polynomial  Matrices 


This  algorithm  is  similar  to  the  one  above  except  it  performs  the  zeroing  process  in  a 
leading  principal  minor  oriented  fashion  so  that  the  algorithm  consists  of  N  —  1  stages 
where  the  k  x  k  leading  principal  submatrix  is  in  a  triangular  form  by  the  end  of  the 
kth  stage.  Furthermore,  the  algorithm  employs  an  additional  substage  which  reduces  the 
degrees  of  the  polynomial  entries  above  the  diagonal  on  the  fly  using  pseudo-division  as  in 
Algorithm  M.  The  order  in  which  the  degrees  are  reduced  is  important  and  is  based  upon 
notions  from  [16]  for  triangularizing  matrices  in  M[Z],  The  order  is  shown  pictorially 
below. 
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The  output  matrix  is  in  column  integral- Hermite  form,  not  simply  triangularized  as  in 
Algorithm  T ,  but  with  the  entries  above  the  diagonal  of  degree  less  than  the  diagonal  entry. 
Clearly,  the  monic  column  Hermite  form  is  easily  obtained  by  left  multiplication  with  the 
appropriate  diagonal  matrix  of  rational  numbers,  a  unimodular  matrix  with  respect  to 
M[Q[s)}. 


BEGIN: 

For  k  =  2  thru  N  Do 

For  n  =  1  thru  k  —  1  Do  (triangularize  the  k  x  kih  leading  principal  minor) 

If  deg  An,n  >  deg  Ak,n  Then  Ak,.  <-»  An>. 

EndlessLoop 

num  <—  pseudoquotient(Ak:m  An  n) 
denom  <—  pseudoremainder(Ak,n ,  An,n) 

Ak,.  <—  denom  *  Ak,.  —  num  *  An,. 

Ak,.  <—  Ak,J GC D{content(Ak,i),  ■ .  ■  ,content(Ak ,  iv)} 

If  Ak,n  7^  0  then  Ak,.  <->  An else  Exit  EndlessLoop 
End  EndlessLoop 
End  Do 

For  i  —  —  1  thru  —k  +  1  step  —1  Do  (reduce  degs  of  abv  diag  polys  in  k  x  kih  minor) 
For  j —  i  +  1  thru  0  Do 

If  degNt+!ifc+J  >  degAk+j,k+j  Then 

num  <- pseudoquotient(Ak+i,k+j,Akk+j,kk+j) 
denom  <—  pseudoremainder (Ak+i,k+i ,  Akk+j,kk+j) 

Ak+i,.  denom  *  .4  —  num  *  Ak+j,. 

Ak+i <—  Ak+i,./GCD{coritent(Ak+i,i), ,  contcnt(Ak+i,N )} 


19 


Endlf 


EndDo 

EndDo 

EndDo 

End 

We  close  this  section  by  noting  that  in  both  Algorithm  T  and  Algorithm  P  each 
pseudo-Euclidean  division  step  affects  the  entire  row  and  the  row  content  is  removed  after 
each  division  step  .  Alternatively,  one  could  solve  a  scalar  Bezout  identity  for  each  zero 
to  be  introduced  using  pseudo-division  techniques  and  then  perform  a  single  elementary 
row  operation  followed  by  a  the  row  content  is  removedrow  content  removal.  However,  the 
single  row  content  of  the  latter  method  will  be  much  larger  than  any  of  the  “elementary” 
row  contents  computed  by  Algorithm  T  or  Algorithm  P.  This  makes  the  alternative  method 
much  less  attractive  than  at  first  glance  in  light  of  the  fact  that  computing  the  many  “small” 
row  contents  is  more  efficient  than  computing  the  single  “large”  row  content. 

8.  Simulation  Results 

Simulations  were  performed  to  determine  the  average  time  required  to  triangularize 
a  square  polynomial  matrix  and  the  maximum  coefficient  length  of  the  output  matrix 
using  both  Algorithm  T  and  Algorithm  P  (see  attached  graphs).  The  maximum  coefficient 
length  is  the  number  of  digits  of  the  largest  (in  absolute  value)  coefficient  appearing  in 
any  polynomial  entry  of  the  output  matrix.  Each  matrix  had  polynomial  entries  with 
randomly  generated  integer  coefficients  uniform  on  [—99,  99].  Runs  were  parameterized  by 
the  dimension  of  the  matrix,  which  ranged  from  2  to  16  and  the  maximum  degree  of  its 
polynomial  entries,  chosen  uniformly  on  [0,  degreemax],  as  degreemax  ranged  from  1  to  6. 

These  simulations  were  conducted  on  a  Texas  Intruments  Explorer  II  with  16  mb  of 
physical  memory  and  128  mb  of  virtual  memory  running  at  40  MHz  using  the  MACSYMA 
version  of  our  algorithms.  The  graphs  represent  the  results  of  the  simulations  averaged 
over  5  runs.  The  results  indicate  that  Algorithm  T  was  moderately  faster  than  Algorithm 
P  in  triangularizing  matrices  up  to  9  X  9.  At  that  point  -Algorithm  T  was  still  faster  for 
triangularizing  matrices  with  lower  degree  polynomials;  but  slower  in  the  higher  degree 
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polynomials.  This  can  be  attributed  to  the  fact  that  Algorithm  P  requires  less  memory 
during  computations  due  to  its  substage  which  reduces  the  degrees  of  the  polynomials 
above  the  diagonal  on  the  fly.  Therefore  costly  garbage  collections,  a  technique  of  freeing 
dynamically  allocated  memory,  are  reduced. 

It  appears  that  both  of  these  algorithms  run  close  to  exponential  time.  The  slopes 
of  the  semi-log  plots  of  the  timings  increase  slightly  with  increasing  polynomial  degree. 
The  maximum  coefficient  length  was  approximately  the  same  for  each  algorithm  and  the 
coefficient  growth  appears  to  be  sub-exponential  with  increasing  matrix  dimension.  A 
16  x  16  matrix  with  degree  6  polynomials  is  the  largest  that  has  been  attempted  with 
Algorithm  P.  It  required  40  hours  to  triangularize  with  the  resulting  matrix  having  a 
maximum  coefficient  length  of  2115  digits. 

Although  Algorithm,  T  was  faster  than  Algorithm  P  on  the  smaller  matrices,  it  did  not 
have  the  overhead  of  putting  the  matrix  into  a  canonic  form  in  the  process;  Algorithm  P 
transforms  the  input  matrix  into  the  canonic  integral-Hermite  form  as  described  earlier. 
The  output  matrix  of  Algorithm  T  therefore  requires  the  application  of  an  auxiliary  algo¬ 
rithm  to  reduce  the  degree  of  the  polynomial  entries  above  the  diagonal  in  order  to  put  it 
in  strict  integral-Hermite  form.  Of  course  this  is  not  necessary  if  one  is  only  interested  in 
rank  information. 

If  one  keeps  in  mind  the  fact  that  our  simulation  results  were  run  on  full,  random  ma¬ 
trices,  which  tend  to  yield  worst-case  performance,  then  these  simulations  indicate  that  our 
algorithms  in  their  current  state  are  ideally  suited  for  problems  in  which  max{m,n}  <  9. 
Such  problems  include  many  practical  control  system  designs,  textbook  problems  in  a 
classroom/lab  environment,  and  empirical  error  analyses  involved  in  research  for  alterna¬ 
tive  apporaches  to  the  machine  computation  of  triangular  forms  of  polynomial  matrices 
based  on  other  arithmetics  such  as  floating-point  or  residue  arithmetic  [9].  For  larger 
problems,  our  code  can  be  modified  in  various  ways  to  yield  approximate  results  in  much 
less  time  while  providing  some  degree  of  error  control.  For  instance,  after  the  integer  co¬ 
efficients  have  reached  a  certain  prespecified  maximum  size,  the  triangularization  can  be 
interrupted  momentarily  and  the  matrix  A(s)  in  M[Z[s]]  at  its  current  state  of  triangu¬ 
larization  can  be  converted  to  an  associated  matrix  A'(s)  in  M[Q[s]]  by  premultiplication 
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with  a  diagonal  matrix  in  M[Q\.  The  matrix  A/(s)  can  then  be  “floated”  to  any  desired 
decimal  precision  and  then  re-expressed  as  a  matrix  in  M [Q[s]]  and  finally  converted  back 
to  Af  |Z[s]]  to  continue  the  triangularization.  An  ad  hoc  technique  such  as  this  is  certainly 
approximate  but  if  done  properly  can  yield  better  results  than  the  ad  hoc  floating-point 
techniques  currently  used.  Refinements  of  this  idea  for  Algorithm  P  with  error  bounds 
and  simulation  results  will  be  appear  elsewhere. 

We  also  compared  our  Hermite  algorithm  to  the  built-in  Hermite  algorithm  included 
with  the  Scratchpad  II  and  Maple  computer  algebra  packages.  On  a  5  x  5  example  generated 
randomly  as  above  our  code  ran  over  100  times  faster. 
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Time  to  Triangularize  (sec)  -  Column  Oriented  Algorithm 
Polynomial  Degrees  1  thru  6 


o  r 


Square  Matrix  Dimension 


Time  to  Triangul arize  (sec)  -  Minor  Oriented  Algorithm 
Polynomial  Degrees  1  thru  6 
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Maximum  Coefficient  Length  (#  of  digits)  =  Column  Oriented  Algorithm 

Potynomial  Degrees  1  thru  6 


Square  Matrix  Dimension 


Square  Matrix  Dimension 


9.  Summary  of  Functions 


MACSYMA  is  a  Lisp  based  system  for  performing  formal,  symbolic  calculations  using 
both  error-free  and  arbitrary  artihmetic.  Since  it  is  Lisp  based,  MACSYMA  runs  the 
fastest  on  Lisp  machines  -  computers  which  have  Lisp  instructions  hard-coded  into  their 
microprocessors,  such  as  the  Texas  Instruments  Explorer  II.  Mathematica  is  a  relatively  new 
computer  algebra  system  written  in  ‘C’  which  has  capabilities  similar  to  MACSYMA.  We 
have  been  running  Mathematica  on  a  NeXT  machine.  If  an  input  matrix  has  polynomials 
with  decimal  or  floating-point  coefficients,  the  coefficients  must  first  be  re-expressed  as 
ratios  of  coprime  integers — both  MACSYMA  and  Mathematica  have  efficient,  built-in 
functions  to  do  this. 

The  following  is  a  summary  of  the  high-level  auxiliary  programs  which  we  have  to 
date  implemented  in  MACSYMA  and  Mathematica.  They  perform  most  of  the  common, 
high-level  tasks  arising  in  the  frequency-domain  approach  to  control  system  synthesis. 

•  RightM atrix Fractional! (s))  —  Computes  a  right  matrix  fraction  description  of  the 
transfer  function  matrix  H(s),  i.e.,  computes  the  matrices  N(s ),  D(s)  such  that  H (s)  = 
N(s)  D(s)~l.  The  LeftMatrixFraction  description  is  analogously  computed. 

•  Bezout(N(s),  D(s))  —  Finds  the  homogenous  and  particular  solutions  to  the  Be- 

zout  equation,  i.e.,  finds  polynomial  matrices  Xp(s),  Yp(s),  Yh(s)  such  that 

A p(s)  A(s)  +  Yp(s)  N(s)  =  I  and  AT(-s)  X(s)  +  Yh(s )  N(s)  =  0.  Used  for  designing 
feedback  compensators  in  the  frequency  domain. 

•  ColumnReduce(D(s ))  —  Column  reduces  the  polynomial  matrix  D(.s),  i.e.,  multiplies 
D(s)  by  an  appropriate  unimodular  matrix  such  that  the  matrix  of  leading  coefficients 
of  its  entries  is  nonsingular.  RowReduce  is  analogously  computed. 

»  Controller^H^s ))  —  Finds  a  controller  form  realization  of  the  transfer -function  ma¬ 
trix  H(s).  Controllability,  Observer  and  Observability  realizations  are  analogously 
computed. 

«  Hermite(N(s ))  —  Finds  the  canonic  column  Hermite  form  of  the  polynomial  matrix 


•  RightCoprime(N(s),  D(s))  —  Determines  the  greatest  common  right  divisor  of  the 
polynomial  matrices  N(s)  and  D(sf  If  it  is  not  unimodular,  it  is  factored  out  of 
both  matrices  making  them  right  coprime.  Used  for  finding  minimal  realizations. 
Left  Coprime  is  analogously  computed. 

•  Smith(N(s ))  —  Finds  the  Smith  form  of  the  polynomial  matrix  N(s).  This  is  a 
canonic,  diagonal  form  of  a  polynomial  matrix.  SmithM  cMillan(H(s ))  —  Finds  the 

r  Smith-McMillan  form  of  the  transfer  function  matrix  H(s).  This  is  a  canonic,  rational, 
diagonal  form  of  a  matrix  whose  entries  are  ratios  of  polynomials. 
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